

pakkettenlijst <- c("openxlsx","sciplot","car")
lapply(pakkettenlijst,function(x){library(x,character.only = TRUE)})


# LOAD DATA
setwd("")
data <- read.xlsx("")



#### DATA CLEANING #####
# timer3 = page submit

tijd_3<- data$ACM_timer_3 + data$IGJ_timer_3 + data$ANVS_timer_3
mean(tijd_3,na.rm=T)/3

tapply(tijd_3,data$attention2,mean,na.rm=T)

tapply(tijd_3,data$attention2,range,na.rm=T)

# select only those with correct attention checks.
#data <- data[which(data$attention2==1),]


data$trust <- (data$agencytrust_1+data$agencytrust_2+data$agencytrust_3) / 3
data$trust_sector <- (data$sectortrust_1+data$sectortrust_2+data$sectortrust_3+data$sectortrust_4) /4

data$allecondities[data$control==1] <- "Neutrale conditie" 
data$allecondities[data$no_transp==1] <- "Non-transparantie" 
data$allecondities[data$gen_eq_pos==1] <- "Positief frame - generiek" 
data$allecondities[data$gen_eq_neg==1] <- "Negatief frame - generiek" 
data$allecondities[data$gen_episodic==1] <- "Episodisch - generiek" 
data$allecondities[data$gen_bench_social==1] <- "Sociale benchmark - generiek" 
data$allecondities[data$gen_bench_hist==1] <- "Historische benchmark - generiek" 
data$allecondities[data$spec_eq_pos==1] <- "Positief frame - specifiek" 
data$allecondities[data$spec_eq_neg==1] <- "Negatief frame - specifiek" 
data$allecondities[data$spec_episodic==1] <- "Episodisch - specifiek" 
data$allecondities[data$spec_bench_soc==1] <- "Sociale benchmark - specifiek" 
data$allecondities[data$spec_bench_his==1] <- "Historische benchmark - specifiek" 
data$allecondities[data$gen_eq_fine==1] <- "Boete - generiek" 
data$allecondities[data$spec_eq_fine==1] <- "Boete - specifiek" 


#### DESCRIPTIVES #######

# age
mean(data$Leeftijd,na.rm=T)
sd(data$Leeftijd,na.rm=T)

#Dutch population mean age = 42.2
t.test(data$Leeftijd[1:5305],Mu=42.2)

# balance check

oneway.test(data$Leeftijd ~ a)
tapply(data$Leeftijd, a, mean,na.rm=T)



# sex
table(data$sex) # 1 = man, 2 = vrouw
table(a, data$sex[data$spec_eq_fine != 1 & data$gen_eq_fine != 1])

geslacht <- matrix(c(53.99,54.84,52.38),ncol = 3)

chisq.test(geslacht,p=c(1/3,1/3,1/3))


# wel of niet? 

a[data$spec_eq_fine==1] <- "Fine" 
a[data$gen_eq_fine==1] 


# education
# low - mid - high according to statistics Netherlands

# https://www.cbs.nl/nl-nl/nieuws/2019/33/verschil-levensverwachting-hoog-en-laagopgeleid-groeit/opleidingsniveau

# 1+2 = Basisonderwijs Lager/voorbereidend beroepsonderwijs (lbo/vmbo) 
# 3+4+5 = Middelbaar beroepsonderwijs mbo), Hoger algemeen voortgezet onderwijs (havo), Voorbereidend wetenschappelijk onderwijs (vwo)
# 6+7 = Hoger beroepsonderwijs (hbo), Wetenschappelijk onderwijs (wo) 

#Population
#Low = 28.3
#Mid = 38.7
#High = 33.0

table(data$education)

data$education_recoded <- recode(data$education,'1=1;2=1;3=2;4=2;5=2;6=3;7=3')
edu <- table(data$education_recoded)
round(prop.table(edu)*100,2)
x <-c(13.06,44.79,42.15)
chisq.test(x,p=c(28.3,38.7,33.0)/100)

# balance check

oneway.test(data$education ~ a)
tapply(data$education, a, mean,na.rm=T)


# knowledge

know_pretest <- (data$know_pretest_1 + data$know_pretest_2 + data$know_pretest_3 + data$know_pretest_4 + data$know_pretest_5 ) /5
know_posttest <- (data$know_posttest_1 + data$know_posttest_2 + data$know_posttest_3 + data$know_posttest_4 + data$know_posttest_5 ) /5

mean(know_pretest,na.rm = T)
sd(know_pretest,na.rm = T)

mean(know_posttest,na.rm = T)
sd(know_posttest,na.rm = T)

t.test(know_pretest,know_posttest)


#Balance test
oneway.test(know_pretest ~ a)
tapply(know_pretest, a, mean,na.rm=T)

oneway.test(know_posttest ~ a)
tapply(know_posttest, a, mean,na.rm=T)



# politics selfplacme
mean(data$politics_selfplaceme,na.rm = T)
sd(data$politics_selfplaceme,na.rm = T)

oneway.test(data$politics_selfplaceme ~ a)
tapply(data$politics_selfplaceme, a, mean,na.rm=T)



##### H1A en H1B #######



# recode into three groups for H1 and H2
# note + fines are left out and result in NA's
a <- NA
a[data$no_transp==1] <- "Non-transparency cue" 
a[data$control==1] <- "Control condition" 
a[data$gen_eq_pos==1] <- "Transparency" 
a[data$gen_eq_neg==1] <- "Transparency" 
a[data$gen_episodic==1] <- "Transparency" 
a[data$gen_bench_social==1] <- "Transparency" 
a[data$gen_bench_hist==1] <- "Transparency" 
a[data$spec_eq_pos==1] <- "Transparency" 
a[data$spec_eq_neg==1] <- "Transparency" 
a[data$spec_episodic==1] <- "Transparency" 
a[data$spec_bench_soc==1] <- "Transparency" 
a[data$spec_bench_his==1] <- "Transparency" 
# fine = missing
a[data$spec_eq_fine==1] <- "Fine" 
a[data$gen_eq_fine==1] <- "Fine"

# compute ci's instead of sd's
require(lsr)
cfi <- function(x){
  x<- ciMean(x,na.rm = T)
}

# remove missings because of "fine"
trust <- data$trust[a != "Fine"]
a <- a[a != "Fine"]

tapply(trust, a, mean,na.rm=T)
tapply(trust, a, sd,na.rm=T)


# color map
require(ggsci)


# figure treatments
pdf(file = "H1 en H2.pdf",width = 6,height = 6)
par(mfrow=c(1,1), mar=c(5,3,1,1))  #bottom, left,top, right
bargraph.CI(a,trust, # trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            cex.axis=1.5,
            cex.lab=1.5,
            col=pal_locuszoom("default",alpha = 1)(3),
           xlab="Experimental condition",
ylab="Trust"
            )

#bargraph.CI(a,data$trust_sector, # trust in sector
#            las = 1,
#            ci.fun=cfi,
#            ylim=c(3,4),
#            xlab="Experimental condition",
#            ylab="Trust"
#)


dev.off()

### exploratieve analyse per toezichthouder
trust_acm <- data$trust_ACM[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_igj <- data$trust_IGJ[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_anvs <- data$trust_ANVS[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 

figuur <-cbind(c(mean(trust[a=="Control condition"],na.rm=T),mean(trust[a=="Non-transparency cue"],na.rm=T),mean(trust[a=="Transparency"],na.rm=T)), # Overall
      c(mean(trust_acm[a=="Control condition"],na.rm=T),mean(trust_acm[a=="Non-transparency cue"],na.rm=T),mean(trust_acm[a=="Transparency"],na.rm=T)), #ACM
      c(mean(trust_igj[a=="Control condition"],na.rm=T),mean(trust_igj[a=="Non-transparency cue"],na.rm=T),mean(trust_igj[a=="Transparency"],na.rm=T)), #IGJ
      c(mean(trust_anvs[a=="Control condition"],na.rm=T),mean(trust_anvs[a=="Non-transparency cue"],na.rm=T),mean(trust_anvs[a=="Transparency"],na.rm=T)) #ANVS
      )
figuur <- as.matrix(figuur)
colnames(figuur) <- c("Overall","Consumer rights","Health care safety","Nuclear safety")
rownames(figuur) <-  c("Control condition","Non-transparency cue","Transparency")

# draai de matrix (transpose)
figuur <- t(figuur)

barplot(height =figuur, beside = T,col = pal_locuszoom("default",alpha = 1)(4), las=1,ylim=c(3,3.6))

vertrouwen <- c(trust_acm,trust_igj,trust_anvs,trust)
conditie <- rep(c(a),4)
groep <- c(rep("Consumer rights",13882),rep("Health care safety",13882),rep("Nuclear safety",13882),rep("Overall",13882))
groep <- factor(groep,levels = c("Overall","Consumer rights","Health care safety","Nuclear safety"))

pdf(file = "per toezichthouder A.pdf",width = 10,height = 6)
par(mfrow=c(1,1), mar=c(5,5,5,5),xpd=T)
bargraph.CI(groep,vertrouwen,group = conditie,
            las=1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            legend = F,
            col=pal_locuszoom("default",alpha = 1)(3),
            xlab="Experimental condition",
            ylab="Trust in regulatory agency"
)
legend("topright",inset = -.01, c("Control","Non-transparency cue","Transparency"),fill = pal_locuszoom("default",alpha = 1)(3),box.lty=0)
dev.off()

# B
pdf(file = "",width = 10,height = 6)
par(mfrow=c(1,1), mar=c(5,5,5,5),xpd=T)
bargraph.CI(conditie,vertrouwen,group = groep,
            las=1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            #legend = T,
            col=pal_locuszoom("default",alpha = 1)(4),
            xlab="Experimental condition",
            ylab="Trust"
)
legend("topright",inset = -.10, c("Overall","Consumer rights","Health care safety","Nuclear safety"),fill = pal_locuszoom("default",alpha = 1)(4),box.lty=0)
dev.off()



##### H2 ######

trust_sector <- data$trust_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0]
leeftijd <- 2021-data$birth
leergedrag <- know_posttest - know_pretest

# # DEZE
# model1 <- lm(trust ~ a)
# model1 <- lm(data$trust[a !="Fine"] ~ a[a !="Fine"])
# summary(model1)
# 
# durbinWatsonTest(model1)
# 
# stats::confint(model1)
# car::vif(model1) # allen dicht bij 1
# 
# # EN DEZE
# # met controls
# model1.1 <- lm(data$trust[a != "Fine"] ~ a[a != "Fine"] + leeftijd[a != "Fine"]+ as.factor(data$sex[a != "Fine"]) + data$education[a != "Fine"] + leergedrag[a != "Fine"]+ data$politics_selfplaceme[a !="Fine"])
# summary(model1.1)
# plot(model1.1) 
# durbinWatsonTest(model1.1)
# 
# car::vif(model1.1) # allen dicht bij 1
# stats::confint(model1.1)


###### hypothese 2#######
#Citizen trust in the regulatory agency strengthens (moderates) the effect of regulatory transparency and citizen trust in the regulated sector

#regulatory transparency --> trust in sector
#                         ^
#                         |
#                 Trust in agency

trust_sector <- data$trust_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0]
a <- a[data$spec_eq_fine ==0 & data$gen_eq_fine ==0]
a <- factor(a,c("Control condition","Non-transparency cue","Transparency"))
a <- relevel(a,ref="Control condition")

# BASE MODEL ## MODEL I ###
h1model<- lm(trust_sector ~ a)
summary(h1model)

stats::confint(h1model)

# ADD TRUST IN AGENCY ## MODEL II #
h1model1 <- lm(trust_sector ~ a + trust)
summary(h1model1)

# Assumptions
durbinWatsonTest(h1model1)
vif(h1model1)

stats::confint(h1model1)

# interaction model OVERALL # MODEL III ##

# centering

centered_trust <- trust - mean(trust,na.rm=T)
conditieXtrust <- a * centered_trust

h1model3 <- lm(trust_sector ~ a * centered_trust)
summary(h1model3)
stats::confint(h1model3)
# Assumptions
durbinWatsonTest(h1model3)
vif(h1model3,type ="high-order") # voor de interactie term

# interacties per sector organisatie
# DV/vertrouwen in sector
trust_acm_sector <- data$trust_ACM_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_igj_sector <- data$trust_IGJ_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_anvs_sector <- data$trust_ANVS_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
# IV / vertrouwen in toezichthouder
trust_acm <- data$trust_ACM[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_igj <- data$trust_IGJ[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_anvs <- data$trust_ANVS[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 


# ACM 
# vertrouwen in toezichthouder
centered_trust_ACM <- trust_acm - mean(trust_acm,na.rm=T)

h1model3_ACM <- lm(trust_acm_sector ~ a * centered_trust_ACM)
summary(h1model3_ACM)
stats::confint(h1model3_ACM)


# IGJ
centered_trust_IGJ <- trust_igj - mean(trust_igj,na.rm=T)

h1model3_IGJ <- lm(trust_igj_sector ~ a * centered_trust_IGJ)
summary(h1model3_IGJ)
stats::confint(h1model3_IGJ)

# ANVS

centered_trust_ANVS <- trust_anvs - mean(trust_anvs,na.rm=T)

h1model3_ANVS <- lm(trust_anvs_sector ~ a * centered_trust_ANVS)
summary(h1model3_ANVS)
stats::confint(h1model3_ANVS)




# controle variabelen

leeftijd[a != "Fine"]+ as.factor(data$sex[a != "Fine"]) + data$education[a != "Fine"] + leergedrag[a != "Fine"]+ data$politics_selfplaceme[a !="Fine"]


h1model2 <- lm(trust_sector ~ a + trust + leeftijd[data$spec_eq_fine ==0 & data$gen_eq_fine ==0]+
                 data$sex[data$spec_eq_fine ==0 & data$gen_eq_fine ==0]+
                 data$education[data$spec_eq_fine ==0 & data$gen_eq_fine ==0] + 
                 data$politics_selfplaceme[data$spec_eq_fine ==0 & data$gen_eq_fine ==0])
summary(h1model2)

# Assumptions
durbinWatsonTest(h1model2)
vif(h1model2,type ="high-order") # voor de interactie term





#oud

require(car)

#leveneTest(trust~a) # significance = not equal. Now equal. 
#H1 <- aov(trust~a)
#summary(H1)
# bonferroni
#(0.05/6) > summary(H1)[[1]][[1,"Pr(>F)"]]
#TukeyHSD(H1)


# effect size
# eta^squared = treatment SSR / SSR Total
# require(sjstats)
options(scipen=999)
# effectsize::eta_squared(H1)
# round(8.73e-03)
# round(0.00873,2)







#### H3 ####
# negative vs positive
h3 <- NA

h3[data$gen_eq_pos==1] <- "Positive framing" 
h3[data$gen_eq_neg==1] <- "Negative framing" 
h3[data$spec_eq_pos==1] <-"Positive framing" 
h3[data$spec_eq_neg==1] <- "Negative framing" 

pal1 <- c("#5cB85CFF","#46B8DAFF")

# negative vs positive



pdf(file = "H2 - H5.pdf",width = 10,height = 10)

par(mfrow=c(2,2),mar=c(6,6,4,6)) #bottom, left, top, right
bargraph.CI(h3,data$trust_sector, # trust in agency
            las = 1,
           ci.fun=cfi,
            ylim=c(3,3.5),
            xlab="Experimental condition",
            ylab="Trust in regulated sector",
          col=pal2,
          main = "H2: Targeted transparency that includes a negative equivalence frame \n has a negative effect on citizen trust in the regulated sector, \n compared with a positive equivalence frame",
            cex.main = 1
)


# 4) Negative episodic information, compared with statistical information from the positive and negative equivalence frame, has a negative effect on citizen trust in the regulated sector. 
h4 <- NA
h4[data$gen_eq_neg==1] <- "Negative statistical" 
h4[data$gen_episodic==1] <- "Negative anecdotal" 
h4[data$spec_eq_neg==1] <- "Negative statistical" 
h4[data$spec_episodic==1] <- "Negative anecdotal" 

pal2 <- c("#357EBDFF","#EEA236FF")

bargraph.CI(h4,data$trust_sector, # trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.5),
            xlab="Experimental condition",
            ylab="Trust in regulated sector",
           col= pal2,
            main = "H3: Targeted transparency that includes negative anecdotal information \n has a negative effect on citizen trust in the regulated sector, \n compared with statistical information.",
            cex.main = 1
)


# H5 social vs historical

h5 <- NA

h5[data$gen_bench_social==1] <- "Social" 
h5[data$gen_bench_hist==1] <- "Historical" 
h5[data$spec_bench_soc==1] <- "Social"
h5[data$spec_bench_his==1] <- "Historical"

bargraph.CI(h5,data$trust_sector, # trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.5),
            xlab="Experimental condition",
            ylab="Trust in regulated sector",
           col= pal2,
            main = "H4: Targeted transparency that includes social reference points \n have a negative impact on citizen trust in the regulated sector, \n compared with historical reference points.",
            cex.main = 1
)



# generic vs specific

h6 <- NA

h6[data$gen_eq_pos==1] <- "Abstract"
h6[data$gen_eq_neg==1] <- "Abstract" 
h6[data$gen_episodic==1] <- "Abstract"
h6[data$gen_bench_social==1] <- "Abstract"
h6[data$gen_bench_hist==1] <- "Abstract"

h6[data$spec_eq_pos==1] <- "Specific"
h6[data$spec_eq_neg==1] <- "Specific"
h6[data$spec_episodic==1] <- "Specific"
h6[data$spec_bench_soc==1] <- "Specific"
h6[data$spec_bench_his==1] <- "Specific"

bargraph.CI(h6,data$trust_sector, # trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.5),
            xlab="Experimental condition",
            ylab="Trust in regulated sector",
          col=pal2,
            main = "H5: Targeted transparency that includes specific information \n has a negative effect on citizen trust in the regulated sector, \n compared with abstract regulatory transparency.",
            cex.main = 1
)


dev.off()



# H3
tapply(data$trust_sector, h3, mean, na.rm=T)
tapply(data$trust_sector, h3, sd, na.rm=T)

b <- t.test(data$trust_sector~h3)
b
require(lsr)
effsize::cohen.d(data$trust_sector ~ h3)

# bonferroni
(0.05/6) > b$p.value




#### H4 ####


tapply(data$trust_sector, h4, mean, na.rm=T)
tapply(data$trust_sector, h4, sd, na.rm=T)

c <- t.test(data$trust_sector~h4)
c
require(lsr)
effsize::cohen.d(data$trust_sector ~ h4)

# bonferroni
(0.05/6) > c$p.value


# H5
tapply(data$trust_sector, h5, mean, na.rm=T)
tapply(data$trust_sector, h5, sd, na.rm=T)

d <- t.test(data$trust_sector~h5)
d

# H6

tapply(data$trust_sector, h6, mean, na.rm=T)
tapply(data$trust_sector, h6, sd, na.rm=T)

c <- t.test(data$trust_sector~h6)
c
require(lsr)
effsize::cohen.d(data$trust_sector ~ h6)

# bonferroni
(0.05/6) > c$p.value

#### FIGUUR MET VERSCHILLENDE TOEZICHTHOUDERS#####

# voor de verschillende sectoren
trust_sector_acm <- data$trust_ACM_sector#[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_sector_igj <- data$trust_IGJ_sector#[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ] 
trust_sector_anvs <- data$trust_ANVS_sector#[data$spec_eq_fine ==0 & data$gen_eq_fine ==0 ]

h3 <- NA

h3[data$gen_eq_pos==1] <- "Positive framing" 
h3[data$gen_eq_neg==1] <- "Negative framing" 
h3[data$spec_eq_pos==1] <-"Positive framing" 
h3[data$spec_eq_neg==1] <- "Negative framing" 



vertrouwen_sector <- c(trust_sector_acm[!is.na(h3)],trust_sector_igj[!is.na(h3)],trust_sector_anvs[!is.na(h3)],data$trust_sector[!is.na(h3)])
conditie1 <- c(rep(c(h3[!is.na(h3)]),4),rep())
groep1 <- c(rep("Consumer rights",4076),rep("Health care safety",4076),rep("Nuclear safety",4076),rep("Overall",4076))
#groep1 <- factor(groep,levels = c("Overall","Consumer rights","Health care safety","Nuclear safety"))

by(vertrouwen_sector,list(groep1,conditie1),FUN = mean,na.rm=T)

pdf(file = "per sector FIGUUR 5.pdf",width = 14,height = 10)

par(mfrow=c(2,2), mar=c(5,5,5,5),xpd=T, cex.lab=1.2)
asgrootte <- 1.2 # textgrootte assen
titelgrootte <- 1.2 # textgrootte
legendcex <- 1.2
labelgrootte <- 1.2
xlabelgrootte <- 1.1
## 1
bargraph.CI(groep1,vertrouwen_sector,group = conditie1,
            las=1,
            ci.fun=cfi,
            #function(x) {c(mean(x)-1.96*se(x), mean(x)+1.96*se(x))}  # 90% CI's
            ylim=c(3,3.6),
            
            #legend = T,
            col=pal2,
            xlab="Experimental condition",
            ylab="Trust in the regulated sector",
            main = "H3: A negative equivalence frame, compared with a positive equivalence \n frame, has a negative effect on citizen trust in the regulated sector.",
            cex.main = titelgrootte,
            cex.axis=asgrootte,
            cex.lab = labelgrootte,
            cex.names = xlabelgrootte
            
)
legend("topright",inset = -.10, c("Negative","Positive"),fill = pal2,box.lty=0,cex=legendcex)

## 2 

h4 <- NA
h4[data$gen_eq_neg==1] <- "Negative statistical" 
h4[data$gen_episodic==1] <- "Negative episodic" 
h4[data$spec_eq_neg==1] <- "Negative statistical" 
h4[data$spec_episodic==1] <- "Negative episodic" 

pal2 <- c("#357EBDFF","#EEA236FF")

# iets andere gemiddelden doordat er niet gedeeld wordt door NA
# opsomming per categorie
by(vertrouwen_sector2,list(groep2,conditie2),FUN = mean,na.rm=T)


vertrouwen_sector2 <- c(trust_sector_acm[!is.na(h4)],trust_sector_igj[!is.na(h4)],trust_sector_anvs[!is.na(h4)],data$trust_sector[!is.na(h4)])
conditie2 <- c(rep(c(h4[!is.na(h4)]),4),rep())
groep2 <- c(rep("Consumer rights",4081),rep("Health care safety",4081),rep("Nuclear safety",4081),rep("Overall",4081))


bargraph.CI(groep2,vertrouwen_sector2, group = conditie2,# trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            xlab="Experimental condition",
            ylab="Trust in the regulated sector",
            col= pal2,
            main = "H4: Negative episodic information, compared with statistical information \n has a negative effect on citizen trust in the regulated sector. ",
            cex.main = titelgrootte,
            cex.axis=asgrootte,
            cex.lab = labelgrootte,
            cex.names = xlabelgrootte
)
legend("topright",inset = -.10, c("Anecdotal","Statistical"),fill = pal2,box.lty=0,cex=legendcex)



## 3
h5 <- NA

h5[data$gen_bench_social==1] <- "Social" 
h5[data$gen_bench_hist==1] <- "Historical" 
h5[data$spec_bench_soc==1] <- "Social"
h5[data$spec_bench_his==1] <- "Historical"

vertrouwen_sector3 <- c(trust_sector_acm[!is.na(h5)],trust_sector_igj[!is.na(h5)],trust_sector_anvs[!is.na(h5)],data$trust_sector[!is.na(h5)])
conditie3 <- c(rep(c(h5[!is.na(h5)]),4),rep())
groep3 <- c(rep("Consumer rights",4083),rep("Health care safety",4083),rep("Nuclear safety",4083),rep("Overall",4083))

by(vertrouwen_sector3,list(groep3,conditie3),FUN = mean,na.rm=T)

bargraph.CI(groep3,vertrouwen_sector3,group = conditie3,
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            xlab="Experimental condition",
            ylab="Trust in the regulated sector",
            col= pal2,
            main = "H5: Negative information about the regulated sector, will have greater negative impact \n on citizen trust when it is framed as social, as opposed to historical, \n comparative information.",
            cex.main = titelgrootte,   
            cex.axis=asgrootte,
            cex.lab = labelgrootte,
            cex.names = xlabelgrootte
)
legend("topright",inset = -.10, c("Historical","Social"),fill = pal2,box.lty=0,cex=legendcex)


## 4
h6 <- NA

h6[data$gen_eq_pos==1] <- "Abstract"
h6[data$gen_eq_neg==1] <- "Abstract" 
h6[data$gen_episodic==1] <- "Abstract"
h6[data$gen_bench_social==1] <- "Abstract"
h6[data$gen_bench_hist==1] <- "Abstract"

h6[data$spec_eq_pos==1] <- "Specific"
h6[data$spec_eq_neg==1] <- "Specific"
h6[data$spec_episodic==1] <- "Specific"
h6[data$spec_bench_soc==1] <- "Specific"
h6[data$spec_bench_his==1] <- "Specific"

vertrouwen_sector4 <- c(trust_sector_acm[!is.na(h6)],trust_sector_igj[!is.na(h6)],trust_sector_anvs[!is.na(h6)],data$trust_sector[!is.na(h6)])
conditie4 <- c(rep(c(h6[!is.na(h6)]),4),rep())
groep4 <- c(rep("Consumer rights",10206),rep("Health care safety",10206),rep("Nuclear safety",10206),rep("Overall",10206))

by(vertrouwen_sector4,list(groep4,conditie4),FUN = mean,na.rm=T)


bargraph.CI(groep4,vertrouwen_sector4,group = conditie4, # trust in agency
            las = 1,
            ci.fun=cfi,
            ylim=c(3,3.6),
            xlab="Experimental condition",
            ylab= "Trust in the regulated sector",
            col=pal2,
            main = "H6: Specific regulatory transparency has a stronger (negative) effect on citizen \n trust in the regulated sector, as opposed to abstract regulatory transparency.",
            cex.main = titelgrootte,
            cex.axis=asgrootte,
            cex.lab = labelgrootte,
            cex.names = xlabelgrootte
)
legend("topright",inset = -.10, c("Abstract","Specific"),fill = pal2,box.lty=0,cex=legendcex)



dev.off()


####### MEDIATION / MODERATION MODEL JPART #####

# http://www.regorz-statistik.de/en/mediation_process_for_r.html#download


# eerst handmatig de process macro inladen!!!

manipulatie <- NA
manipulatie[a=="Transparency"] <- 1 
manipulatie[a=="Control condition"] <- 0
manipulatie[a=="Non-transparency cue"] <- 0
data$manipulatie <- manipulatie

# let op = boete zit erbij. op de plek waar 'a' gemaakt wordt.
# a inladen vanaf regel ~ 160. Niet de data$a maar losse a 

# model 1, zonder boete. Controle & non-transparantie bij elkaar
#transparantie --> vertrouwen sector
#                ^
#                |
#            vertrouwen
#            toezichthouder
            

modeldata <- data.frame(trust_sector = data$trust_sector[data$spec_eq_fine ==0 & data$gen_eq_fine ==0],
                        trust = data$trust[data$spec_eq_fine ==0 & data$gen_eq_fine ==0],
                        manipulatie = data$manipulatie[data$spec_eq_fine ==0 & data$gen_eq_fine ==0])

process(data = modeldata,y="trust_sector",x="manipulatie",w="trust",model = 1)


# model 4, zonder boete. Controle & non-transparantie bij elkaar
#transparantie    -->     vertrouwen sector
#       \                 /  
#        \               /
#          \ vertrouwen /
#            toezichthouder

interactie <- modeldata$trust*modeldata$manipulatie
modeldata[,4] <- interactie
process(data = modeldata,y="trust_sector",x="manipulatie",m="trust",model = 4)





require(mediation)
set.seed(2014)
data("framing",package = "mediation")
# example
?framing
med.fit <- lm(emo~treat + age + educ + gender + income, data = framing)
out.fit <- glm(cong_mesg ~ emo + treat + age + educ + gender + income , data= framing, family = binomial("probit"))

met.out <- mediate(med.fit,out.fit,treat="treat",mediator="emo",robustSE=TRUE, sims=100)
summary(met.out)
plot(met.out)


# For HHG

interactie <- modeldata$trust*modeldata$manipulatie
modeldata[,4] <- interactie
process(data = modeldata,y="trust_sector",x="manipulatie",m="trust",model = 4)


med.fit <- lm(modeldata$trust_sector ~ modeldata$manipulatie)
out.fit <- lm(modeldata$trust_sector ~ modeldata$trust + modeldata$manipulatie)

# ACME = average causal mediation effect
# ADE = average direct effect

met.out <- mediate(med.fit,out.fit,treat="modeldata$manipulatie",mediator="modeldata$trust", robustSE = TRUE, sims=100)
summary(met.out)
plot(met.out)

sed.est <- mediate.sed("trust_sector","trust","manipulatie", SI = TRUE, boot = TRUE, sims=1000,data= modeldata)
summary(sed.est)



